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For large scale electronic structure calculation, the Krylov subspace method is introduced 
to calculate the one-body density matrix instead of the eigenstates of given Hamiltonian. This 
method provides an efficient way to extract the essential character of the Hamiltonian within a 
limited number of basis set. Its validation is confirmed by the convergence property of the den- 
sity matrix within the subspace. The following quantities are calculated; energy, force, density 
of states, and energy spectrum. Molecular dynamics simulation of Si(OOl) surface reconstruc- 
tion is examined as an example, and the results reproduce the mechanism of asymmetric surface 
dimer. 
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1. Introduction 

The study of nanoscale systems requires a large-scale 
atomistic simulations with quantum mechanical free- 
doms of electrons. The practical requirement to carry out 
the simulations is how to extract desired quantities from 
a given large Hamiltonian matrix, not only accurately 
but also efficiently. Simulation methods in large scale 
systems have been studied already in the last decade. -^"'^ 
In order to execute molecular dynamics simulation, one 
needs information about the total energy and forces on 
an individual atom, and these physical quantities should 
be obtained by means of either eigen states | <^q) or the 
one-body density matrix p of the system; 



a 



(1) 



) is the Fermi-Dirac distribution function 



Here 

as a function of the eigen energy of the eigen states 
I (pa) and the chemical potential /i of the system. 

The molecular dynamics calculation in large-scale sys- 
tems can be done on the basis of transferable short-range 
tight-binding Hamiltonians H , where we calculate the 
physical property {X) as 



{X)=TT[pX]=Y,p,,X,,. 



(2) 



Here i and j are suffices of atomic site and orbital. The 
energy and forces acting on an atom are contributed only 
by elements that have non-zero values of the Hamiltonian 
matrix. In other words, even though the density matrix is 
of long range, only the short range behavior of the density 
matrix is essential.'' Therefore, the essential methodology 
for large scale calculations is how to obtain the short 
range part of the density matrix p without calculating 
eigen states of the original Hamiltonian. The essential 
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point here is the fact that we adopt the short-range tight- 
binding Hamiltonian and this makes computation local. 
We will only comment here that the short-range tight- 
binding Hamiltonian can be always constructed both in 
insulators and metals from the first principle theory.^' ^ 

We have developed a set of methods, without calcu- 
lating eigenstates, of large-scale atomistic simulations, 
which are based on generalized Wannier state and hy- 
brid scheme within fully quantum mechanical descrip- 
tion of electron systems. ^"'^ These methods are rigor- 
ously a linear scale simulation in atom number, and 
were tested upto 10^ atoms by using a standard work- 
station. The generalized Wannier states are defined for- 
mally as unitary transformation of the occupied eigen 
states, though eigen states are not actually obtained. 
This method is practical and efficient in covalent bonded 
materials, where the localized Wannier states reproduce 
well the electronic structure energy and the density ma- 
trix, at least its short-range behavior. We observed that 
the bond forming and breaking processes are well de- 
scribed in the localized Wannier states as changes be- 
tween a bonding and non-bonding orbital.^' In metal- 
lic systems, however, situations are quite different and 
other practical methods should be developed. 

The aim of the present work is to establish an novel 
extension of methodology practical in metals. We will 
develop a novel method based on the Krylov subspace 
(KS) method to achieve computational efficiency. In §2 
we review the KS method and the density matrix is rep- 
resented in the KS. An example will be presented based 
on our numerical results. These include a discussion of 
locality of off-diagonal elements of the density matrix. In 
§3, as an example of the molecular dynamics simulation, 
the reconstruction of Si (001) surface will be discussed. 
We will also show how the energy spectrum can be ob- 
tained in our developed method. In §4 we summarize the 
work presented in this paper. 
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2. Density Matrix Calculation based on Krylov 
Subspace Method 

In this section, we will show theoretical background 
of the KS (Krylov subspace) method to extract density 
matrix for molecular dynamics simulation. Short review 
of the KS method is followed by analysis of its arithmetic 
structure including convergence property which justifies 
the present method. 

2.1 Krylov subspace method^^' 

The KS (Krylov subspace) method gives the math- 
ematical foundation of many numerical iterative algo- 
rithms such as the conjugate gradient method. This 
method provides an efficient way to extract the essential 
character of the original Hamiltonian within a limited 
number of basis set. Starting from a certain vector | i), 
a subspace of the original Hilbert space is generated by 
a set of vectors; 

M), H\^).^ H^\^).^ H-^--'\z). (3) 

The subspace spanned by the basis vectors {ff" | i)} in 
eq. (3) is generally called the Krylov subspace (KS) in 
the mathematical textbooks. The dimension of the KS is 
denoted as ^k- We will denote the ortho normalized basis 
vectors in the KS as 

|if«)(^|z)), (4) 

Since the matrix H is Hermitian, the Gram-Schmidt or- 

thonormalization procedure gives one possible (but not 
necessary) choice of the basis set that satisfies the three- 
term recurrence relation called the Lanczos process; 

bn I kI:1,) = {H- an) I - b:_, I if f i), (5) 

with b-i = 0. Hereafter we restrict ourselves to real sym- 
metric Hamiltonian matrix, H. 

From the practical point of view of calculations, the 
procedure of matrix- vector multiplication, H \ Kn^), 
(;onsumes the CPU time mostly, then the number of bases 
in the KS (i/k) should be chosen to be much smaller than 
that of the original Hamiltonian matrix. This drastic re- 
duction of the matrix size or the dimension of the KS is 
a great advantage for a practical large-scale calculations. 
The dimension of the KS i^k should be chosen, for exam- 
ple, as fK = 30. We then denote the reduced Hamiltonian 
as for the KS {| JsT^'^)}- 

2.2 Density matrix calculation in the Krylov subspace 
In order to extract desired density matrix, we diago- 

nalize the reduced Hamiltonian matrix . Once one 

obtains the eigenvalue e^a and eigenvector | w'^) as 

H^ii) I yji^) = 4') I ^i,W), (6) 
the eigen vector can be expanded in terms of the basis 

l^i^^> = J2^^n\K^)- (7) 



Author Name 

We introduce the density matrix operator within the KS: 

The essence of the present method is the replacement 
of the density matrix (« | p | j) by that of the KS {i \ 

I /5 h') ^ I /5^« (9) 

Once this procedure is allowed, it is a great advantage 
from the view point of practical calculations. 
Let us introduce the projection operator; 

pKii) ^ 5^|t„W)(t„«|=5^|/f«)(/f«|(10) 

a n 

The crucial point is for the calculation of (i | /5 | j) that, 
though the state | i) is an element of the KS (P^(') | 
i) =1 i)), the state | )) may be not an element completely 
included within the KS (P^(*) \ j) ^| j) nor 0). Even so, 
the density matrix of the KS {i \ p^(') | j) holds the 
following relation; 

(^ I I j) = ^(^ I I if«)(ifW I j). (11) 

n 

To show eq. (11) we use a relation 

The replacement eq. (9) is rigorous when z^k is equal to 
the dimension of the original Hilbert space. When i^K is 
much smaller, this replacement (9) can be justified only 
if the convergence of the summation in eq. (11) is fast 
enough and the contribution from large n is negligible. 
We will check the n dependence of both {i \ p^^'^ | Kn^) 
and {Kn^ \ j) in Subsection 2.3. 

Considering the spin degeneracy, the relation of elec- 
tron number Neiec and the chemical potential fi can be 
given as 

i 

= (13) 

which is used to determine the chemical potential p in 
the system. 

For short summary of this subsection, we note that 
the essential procedure is only the part reducing the di- 
mension of the original Hamiltonian matrix H to that of 
the reduced matrix H^^^K Once we obtain {| Wa"^)} the 
cost of calculating {{i \ p^'*' | j)} for necessary enough 
number of neighboring sites and orbitals j of a fixed i 
is of the order of one, independent of the system size or 
the total number of atoms. And, furthermore, the calcu- 
lation of them is perfectly parallelizable with respect to 
sites and orbitals i. 

2.3 Convergence properties of the density matrix 

In order to demonstrate the validity of the replacement 
eq. (9), we check the convergence of eq. (11). The con- 
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Fig. 1. Decay properties of the ofT-diagonal density matrix, p^''\ 
of 262,144 atom for Si (solid square with Une) and that for C 
(open square with Une) as a function of n in the summation of 
eq. (11) representing number of hoppings from center atom. The 
dashed Une (two dot dashed line) is a guide to the eye repre- 
senting 1/n (1/ra^) behavior. Inset shows pf^''' for Si in linear 
scale. 

vergence varies according to the locality of the original 
Hamiltonian as well as the choice of starting basis or- 
bitals. To demonstrate the property, we choose a crystal 
of diamond structm'e with 262,144 atoms using a trans- 
ferable tight-binding Hamiltonian.^'* We also choose the 
four sp'^ orbitals per atom as starting basis orbitals 

Figure 1 shows decay property of pf^^\= {i \ p^^'-'i | 
Kn'^)) as a function of the n in the summation of eq. (11). 
As shown in the inset, pfj^^^ decays oscillatery. We plot 
its absolute value to see the decay behavior. In case of 
silicon, we can read that p^^*^ decays as fast as l/n?. 
On the other hand the value of {Kn'^ \ j) also decays 
as a function of n due to the fact that the state | Kn^) 
extends over sites reached with n-steps from the starting 
state I iff^)(=| i)) to cover another localized basis |j) on 
one site. Decay property of {Kn'' \ j) depends on j but its 
maximum value in the present system decays as 1 /n (not 
shown in the figure). Therefore, the products in eq. (11) 
decays as We examined several cases in different 

system size (512, 4096, 32768, and 262144 atoms), and 
found that the decay property is almost independent of 
the system size. In case of carbon, on the other hand, 
the decay rate of p^^'"* is even more faster, which can be 
understood from the locality of the Wannier state. 

Since the choice of the starting basis is arbitrary, we 
can choose the four atomic orbitals at each atom site, (s, 
Pk, Pj/, Pz) , as starting basis orbitals . Here, however, 
we choose the starting bases \i) =\ k'^'') as the four sp'^ 
orbitals, because the cohesive mechanism is clarified with 
such hybridized bases. Due to the crystalline symmetry 
of diamond structure, the four sp'^ bases are equivalent 
and only one example is enough for the explanation of 
the cohesive mechanism. The dominant interaction in the 
Hamiltonian is the hopping along the sp^ bond. If we 
ignore other hoppings in the Hamiltonian, the Krylov 
subspace with Vk = 2 gives the sp'^ bonding and anti- 
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bonding orbitals as 

I K^)± I K^) 

which forms a desirable basis set in the present case. 

We would consider an example of possible slowest con- 
vergence of eq. (11) where we can define the Fermi wave 
vector k-p] Since the three-term recurrence relation in 
eq. (5) suggests a mapping of the original system to a 
one-dimensional chain model, it is instructive to com- 
pare with simple consideration of one dimensional sys- 
tem with constant energy a and hopping b in eq. (5). 
This case corresponds to the one-dimensional free space, 
in the continuum limit, and the density matrix is given 
by analytic form 

pix.x') = r e"=(^-^')d/fc = sinfcF(x-xO ^ ^^g^ 

This can be understood as 1/n behaviour of pi„ with 
oscillation. Even in this case {Kn^ \ j) decays as 1/n 
and the products in eq. (11) decays as 1/n^. Though 
the analysis for other realistic systems like simple cubic 
lattice will be shown elsewhere, we should mention that 
there are several practical examples where p^/*'' decays 
as 1/n. 

Further, from the view point of practical calculations, 
the decay rate can be controlled by the temperature fac- 
tor k-QT; The higher the temperature, the faster the de- 
cay. These facts validate the convergence of the sum- 
mation in eq. (11) and justifies the replacement eq. (9). 



0.11 




Fig. 2. Reduced matrix size dependence of the off-diagonal ele- 
ments of the density matrix, pij, of 512 Si atom (open circle with 
line). The dot dashed Une represents the results of exact diag- 
onalization of the original Hamiltonian. Inset shows schematic 
pictures of the sp'^ hybrid orbitals | hi), | hn), | hm), and | hiv). 
Solid orbitals in each figure represent the combination to con- 
tribute the density matrix. 
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2.4 Convergence properUes of off-diagonal elements of 
the density matrix and the total energy 

As an example of the present method, we show the 
calculated density matrix and compare with that of di- 
agonalization of the original Hamiltonian. We pick up 
two nearest neighbor bond sites along a linear path with 
four sp^ hybrid orbitals | hj), | hn), | hm), | hiv), where 
two orbitals {| hi) and | hn)} and {| hm) and | hjy)} 
are on the same bond sites and {| hn) and | hni)} are 
on the same; atom. Sec inset of Fig. 2 for the configu- 
ration and the phases of respective hybrid orbitals. The 
exact values of these matrix elements are calculated by 
the exact diagonalization of the original Hamiltonian as 
follows; (hi I p I hii) = 0.439, (hn | p \ hm) = 0.078, 
(hi I p I hm) = -0.008, (hi I p I hiy) = -0.071. These 
four are the typical elements between nearest neighbor 
bond sites which can be easily understood from the view 
points of the Wannier states. 

When the dimension of the KS increases, the calcu- 
lated values of off-diagonal elements of the density ma- 
trix gradually approach to the exact values and saturate. 
Figure 2 shows the corresponding results for Si crystal 
with 512 atoms. In the present case they are saturated 
at around i/k = 30. The resultant convergent behavior 
and values are both excellent. 

We note here that the convergence of the to- 
tal energy could not be a unique measure of the 
convergence of the calculations. The convergence of 
the total energy is more rapid in comparison with 
that of the off-diagonal elements of the density ma- 
trix. In fact, the exact value of the band energy is 
—5.082 eV/electron and the calculated deviation from 
this is +80, 4-23, +4, +1, +0 meV/electron for vk = 
7, 10, 20, 25, 30, respectively. 

It must be mentioned that, for the present covalent 
bonded systems, the generalized Wannier state can be 
reasonably reproduced by the first order perturbation 
theory of the sp^ bonding orbitals and the Perturbative 
Order-N method is quite efficient.^' ^ The computational 
cost of the present KS method is less efficient in these 
systems. 

2.5 Computational details and comparison with other 
methods 

In actual computations, we adopt the following proce- 
dure: 

[i] Generate the Krylov subspace defined by eq. (4) and 
generate eigen states within the KS by eq. (6). 

[ii] Determine the chemical potential ji from the diagonal 
elements of the density matrix by using eq. (13). 

[iii] Calculate the off-diagonal elements of the density 
matrix through eqs. (8) and (11). 

[iv] Calculate forces acting on each atom and move 
atoms. 

[v] Return to the procedure [i]. 

The computational time and memory size are mostly 
consumed in the part of generating the Krylov subspace. 
The computational cost of all other procedures is actu- 
ally linearly proportional to the number of atoms. Fur- 
thermore, the only global quantity we use is the chemical 
potential and all other calculation is purely indepen- 



dent with respect to each starting vector. Therefore, the 
computational routine is parallelizable, and actually we 
made use 128 and 256 parallel processors with the Mes- 
sage Passing Interface (MPI) technique. 

Since the present method and the recursion 
method^^"^'' are both based on the construction of 
the Krylov subspace, one might suppose that it were an 
extension of the recursion method. However, it is not the 
case. All calculations in the present method are based 
on the eigen values and eigen vectors in the Krylov 
subspace and one can calculate directly off-diagonal 
elements of the density matrix. On the other hands, the 
recursion method is the way of calculating the diagonal 
Green's function in a form of the continued fraction. 
The discussion in the recursion method is always based 
on the diagonal elements of Green's functions G. The 
proposed way to calculate the off-diagonal Green's 
function in the recursion method may be^^"^*^ 

Gij = 2{pi+j,i+j - Gi-j,i-j^, (16) 

which needs a lot of computational resources. The re- 
cursion method would recommend, in order to calculate 
the off-diagonal Green's function, to use the recurrence 
relation of the Green's function, ^^'^^ but it contains po- 
tential growth of a numerical rounding error. 

The density matrix actually is given by the energy in- 
tegration of the Green's function in the recursion method 
as 

Pij = ~ de ImG,,(£)/(i-^) , (17) 

which causes a numerical error. The present method is 
completely free from above-mentioned difficulties in the 
recursion method. All calculations in the present method 
are based on the eigen values and eigen vectors in the 
Krylov subspace and one can calculate directly diagonal 
and off-diagonal elements of the density matrix simulta- 
neously. 

3. Example : Results and Discussions for the 
Surface Reconstruction of Si (001) 

In this section, we demonstrate how the electronic 
structure within the KS method gives the correct atomic 
structure. We show the results of molecular dynamics 
simulation of Si (001) surface reconstruction of a slab sys- 
tem of 1024 atoms. The essence of the quantum mechani- 
cal freedoms is the fact that sp'^-hybrid bonds are formed 
in the bulk region, but not on surfaces. Specifically sur- 
face atoms move to form asymmetric dimer.^^'^* We will 
show the result of the present method and discuss the 
local electronic structure and the energy spectrum. We 
also examine total energy difference for proposed three 
reconstructed configurations. 

3. 1 Tilt angle of surface dimers 

In ideal Si(OOl) surface, a pair of surface atoms has 
four electrons as dangling bonds. Two of them forms a 
cr-bonding state and a surface dimer appears. The other 
two electrons are directly related to the asymmetric ge- 
ometry of the surface dimer. The Hilbert space for these 
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electrons is restricted to the basis set orthogonal to the 
(T-bonding states and two back-bond states. If the four 
atomic orbitals, (s, Pa;, Py, Pz) per atom are considered, 
three freedoms are excluded by the orthogonality to the 
above three states. In the asymmetric dimer, the re- 
stricted basis set is given by an atomic basis of the upper 
atom with a large s component and a relatively low en- 
ergy level and the one of the lower atom with a large 
p component and a relatively high energy level. Then 
the system can gain the energy, with the increase of s 
component, by charge transfer from the lower atom to 
the upper one. This mechanism is sometimes called 'de- 
hybridization' in the sense that the sp'^-hybridization is 
cancelled (See Ref.^ and the references therein). In our 
previous work, we have observed a dynamical process of 
forming the asymmetric dimer, according to the above 
energy gain mechanism.^ Therefore the present method 
should reproduce the above energy gain mechanism so as 
to reproduce the asymmetric dimer. 

One of the factor to characterize the surface dimer 
is its tilt angle, 6. (See inset of Fig. 3(a).) Theoretical 
and experimental data of tilt angle 9 are reviewed in 
ref.,^^ and are ranging from 5° to 19 ° . The reported 
tilt angle by the exact diagonalization of the same tight- 
binding Hamiltonian is ~ 14 ° while our result based 
on the KS method is 6* = 13.4 ° with the size of the 
KS i^K ~ 30. This result indicates that the present KS 
method extracts the essential character of the original 
Hamiltonian well. We will discuss in the next subsection 
that the asymmetric surface dimer is determined by the 
electronic states close to the chemical potential. 

3.2 Energy spectrum and the density of states 

While methods of density matrix may usually not pro- 
vide an information about energy spectrum of electronic 
structure, the present method can do at the same time. 
To discuss the electronic spectra in the framework of the 
KS method, we introduce the Green's function Gij{e) ; 



G^J{e) = [{e + iS-HY 



(18) 



where S is an infinitcsimally small positive number. Since 
the replacement for the density matrix (9) is guaranteed, 
the same replacement for the Green's function is also 
allowed; 



Gf^<^^(£), 



(19) 



where the matrix elements of the Green's function in the 
KS is defined as 



GaiCa 



e -\-i5 — e, 



(i) 



(20) 



Actually the Green's function Gij{e) can be calculated 
with the Green's function G'^}^\e) in the KS as; 



<«(e) = ^G,^«(e)(i^W|j) 



(21) 



The equation (21) is equivalent to (11) and can be proven 
similarly by using the projection operator P^('). 

In order to single out the physical insight behind the 
asymmetric dimer, we calculate local density of states 
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Fig. 3. {a)Local density of states (IDOS) per atom for the system 
with asymmetric dimer and that for the system of crystal. Solid 
line (broken line) in upper panel represents an upper (lower) 
atom of the asymmetric dimer. (b) COHP and integrated COHP 
for the corresponding dimer. The energy zeroth both in (a) and 
(b) are common and is set to be the top of the occupied states 
in the bulk. In order to show the structure we introduce finite 
imaginary part, 5 = 0.136eV, in the energy denominator of the 
Green function. The size of the reduced matrix is i/k = 30 and 
the temperature factor of the system in eq. (8) is T = 1580 K 
(=0.136eV). The chemical potential is estimated as fi =0.126eV. 



(IDOS) per atom of the system with reconstructed sur- 
face with dimer as shown in Fig. 3(a). The IDOS can be 
defined as 



ni{e) = VlmG,„j„(£) (22) 

TT — ' 

a 

= P<5(£-4'"^), (23) 



a,K 



where / and a are the atomic site and orbitals, respec- 
tively, and K is suffix for eigen states of the KS. First 
of all, we see the IDOS of crystal. Because of the finite 
number of computed levels, i>k = 30, the shown IDOS 
has thirty spikes with weight factor | {la \ Wk"''') ^ 
distributed from bottom to top of the band. Here we 
have introduced finite imaginary part, 6 — 0.136eV 
(10^^ R-yd), to smooth out these spiky structure. The 
calculated IDOS of crystal reproduces the gap that lies 
within ^ leV satisfactory. The IDOS of the deeper 
layer of the present slab system is similar to this and does 
not change before and after the surface reconstruction as 
it should be. In the IDOS for dimerized surface atoms, 
the IDOS of the upper (lower) atom has peak at —1.25 
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(+0.54) eV in Fig. 3(a). The former (latter) peak corre- 
sponds to occupied (unoccupied) surface state and the 
difference of the spectra represents the electron charge 
transfer from the lower atom to the upper atom in the 
asymmetric dimer, as explained in §3.1. In other words, 
the Krylov subspace method reproduces the electronic 
structure in the asymmetric dimcr. 

We note here about the two controlling parameters to 
reproduce the asymmetric dimer; the size of the KS, Vk, 
and the temperature factor of the system, T. Both may 
affect the convergence speed in eq. (11) as well as the 
energy resolution of the simulation. The choice of vk is 
important to reproduce the asymmetric dimer since the 
surface dimer reflects the electronic structure close to the 
chemical potential, in particular the occupied and unoc- 
cupied surface states. The size of the KS should be chosen 
so large that the profile of the surface states are well re- 
produced. Actually, the calculation with vk < 20 leads 
unstable value of 9, for example, 9 = 0.2,9.8. 14.5,4.6 ° 
for z^K = 15, 16, 17, 18, respectively. While those with 
i^K > 25 gives stable value, 13 ~ 14 ° . We have cho- 
sen z^K = 30. The choice of T is also important since 
the surface states are energetically close to the chemical 
potential. The temperature should be chosen so small 
that the occupied and unoccupied surface states are well 
separated energetically. 

In order to see the chemical bonding in condensed mat- 
ters, we introduce the following quantity; 

Cij{e) = -- V lmGic,M^)Hj0j^. (24) 

This is sometimes called the crystal orbital Hamiltonian 
populations (COHP).^'' The integration of this quantity 
gives cohesive energy from a pair of atoms just as the 
integration of local DOS gives occupation number. Ac- 
tually, the total energy is decomposed into contributions 
of each atom pair as a sum of integration over the energy 
oiCij; 

Tv{pH) = J2J2pia,j0Hj0ja (25) 

I, J a,f3 

= E /" (26) 
i,j 

The analysis of the COHP and the integrated COHP 
shows where and how the bond formation stabilizes en- 
ergetically the system. The COHP for the dangling bond 
pair (in ideal surface) is negligible (not zero), because 
interaction matrix element Hj^ja within the dangling 
bond pair is very small due to a larger interatomic dis- 
tance. Once an surface dimer is formed (though the 
atomic pair is the same), the COHP gives a finite value 
(Fig. 3(b)), because the interatomic distance is short- 
ened and the interaction matrix element becomes finite. 
The integration of the COHP has its minimum almost 
at the chemical potential. This is a demonstration of the 
cohesive mechanism of covalent bonded materials. 



3.3 Energy difference between different configurations 
of dimerized Si (001) surface 
The dimers may align on the Si (001) surface with 
three proposed reconstructed surface configurations, (2 x 
1). (2 X 2). and (4 x 2).^'^'^^ Among them, the present 
calculation indicates that the (4 x 2) configuration is 
that of the lowest energy. The calculated energy differ- 
ences from the (4 x 2) structure are 86.7 nieV/dimer 
for -E(2xi) find 0-3 meV/dimer for i?(2x2)- These val- 
ues agree well with the exact calculation using the same 
Hamiltonian, -B(2xi) ~ -£'(4x2) = 73.6 meV/dimer and 
E(2x2) ^ £^(4x2) = 1-2 meV/dinier, respectively.^^ This 
shows that the numerical error with the KS method is 
small and the present method gives a satisfactory results 
in a fine energy scale with tight-binding Hamiltonian. 
On the other hand, we should comment that the tight- 
binding formulation itself can be the another origin of an 
error. In general, the energy scale in meV/atom is too fine 
to discuss in the present tight-binding Hamiltonian. An 
ab initio calculation gives £^(2x1) ~-£'(4x2) = 51± 21 and 
E(2x2) — £-(4x2) = 3± 13 meV/dimer, respectively.^'* 

4. Conclusions 

In the present paper we presented a novel method us- 
ing the Krylov subspace for the molecular dynamics sim- 
ulation based on large-scale electronic structure calcula- 
tion. By means of the reliable treatment of the reduced 
matrix deduced from the Krylov subspace method, the 
method provide an efficient and practical way to calcu- 
late the density matrix. The method also provides a way 
to calculate the energy spectrum on the same standpoint 
as the density matrix. As an example, the method is ap- 
plied to the problem of the surface reconstruction of Si 
(001). We have pointed out through its analysis that the 
appropriate choice of the two controlling parameters, the 
size of the Krylov subspace and the temperature factor, 
is important. Both may affect the computational c;ost 
and the accuracy. Though the present calculation is just 
one example, it leads us to a general guiding principle 
in choosing the optimal values of the controlling param- 
eters. 

In the present methodology the computational proce- 
dure of the density matrix, pij, is independent for each 
atomic orbital, i, except the determination of the chem- 
ical potential, then the present method is very prefer- 
able for the parallel computation. Moreover, this inde- 
pendency of the basis lead us a hybrid scheme within 
quantum mechanics.^ In the hybrid scheme, the density 
matrix is decomposed into sub matrices and the sub ma- 
trices are determined by different methods. Molecular dy- 
namics simulation with 10^ atoms by the hybrid scheme 
between the present KS method and the perturbativc 
Wannier state method is examined and will be published 
elsewhere. 

Since this newly developed method is a general the- 
ory for large matrices, the method is applicable for not 
only covalent bonded materials but also other systems 
like metal. The present KS method has a potentiality of 
wide applicability, even in non-Hermitian matrix, since 
the fundamental concept lies in the general linear algebra 
of large matrices. 
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